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ABSTRACT 

In this paper we examine the issue of characterising the transport associated with 
gravitational instabilities in relatively cold discs, discussing in particular the condi- 
tions under which it can be described within a local, viscous framework. We present 
the results of global, three-dimensional, SPH simulations of self-gravitating accretion 
discs, in which the disc is cooled using a simple parametrisation for the cooling func- 
tion. Our simulations show that the disc settles in a "self-regulated" state, where the 
axisymmetric stability parameter Q k. \ and where transport and energy dissipation 
are dominated by self-gravity. We have computed the gravitational stress tensor and 
compared our results with expectations based on a local theory of transport. We find 
that, as long as the disc mass is smaller than 0.25M+ and the aspect ratio H/R < 0.1, 
transport is determined locally, thus allowing for a viscous treatment of the disc evo- 
lution. 
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; 1 INTRODUCTION 



One of the basic unknowns of accretion disc theory is the 
physical mechanism ultimately responsible for angular mo- 
mentum transport and energy dissipation in the disc. It is 
well known that classical hydrodynamical viscosity is not 
sufficient to provide accretion at the rates inferred from the 
observations in almost every astrophysical context where ac- 
cretion discs play a role. The usual way to overcome this 
difficulty is to assume that transport is dominated by some 
"anomalous" viscous phenomenon, possibly related to col- 
lective instabilities in the disc, and to to give some ad 
hoc parameterisations for the viscosity. The most widely 
used of such parameterisat ion is the so-called a-prescription 
iShakura fc Sunvaevlll973l . see Section 2). 

It has been recently recognised that accretion discs 
threaded by a weak magnetic field ar e subject to MHD insta- 
bilities fsee lBalbus fe HawlevllT99ci and references therein), 
that can induce turbulence in the disc, thereby being able 
to transport angular momentum and to promote the accre- 
tion process. However, in many astrophysically interesting 
cases, such as the outer regions of protostellar discs, the 
ionisation level is expected to be low, reducing significantly 
the effects of magnetic fields in determining the dynamics 
of th e disc, at least over limited radial ranges iGammiel 
1996). A possible alternative source of tra nsport in cold 
discs is provided by gravitational instabilities ijLin fc P ringlc 
ll987tlLaughUn fc Bodenheimerll99llArmitage et al.l200ll) . 



Moreover, in environments where the disc mass is a signifi- 
cant fraction of the central object mass, such as during the 
assembly of protostellar discs, it is inevitable that disc self- 
gravity will play a role. 

The axisymmetric stability of a thin disc with respect 
to its own sel f-gravity is determined by the parameter Q 
jToomreill9'6l) . defined as: 



Q = 



ttGE' 



(1) 



where c s is the effective thermal speed of the disc, E is the 
surface density, and K is the epicyclic frequency (which, for a 
Keplerian disc, is equal to the angular velocity fi). The disc 
is unstable if Q is smaller than a threshold value Q m 1. It 
has been lo ng recognised, especially in the cont ext of galactic 
dynamics fiiohHl97lt iBertin fc Romeolll98Sl) . that the de- 
velopment of gravitational instabilities would lead to a self- 
regulation process: if the disc is initially cold, in the sense 
that Q < Q, then gravitational instabilities would heat it up 
on the fast dynamical timescale, bringing it toward stability; 
on the other hand, if the disc is hot enough to begin with, 
then radiative cooling is going to bring the value of Q down 
toward an unstable configuration. As a result of the pres- 
ence of these two competing mechanisms, the "switch" as- 
sociated with the onset of gravitational instabilities will act 
as a thermostat for the disc, which is therefore expected to 
be always close to marginal stability. A similar approach has 
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been also suggeste d in the case of accretion discs iPaczvriskil 
ll97Sl:lBertirJll997l) . 

From the observational point of view, there are already 
some clues that the disc self-gravity can be important both 
in the context of protostellar discs and in accretion discs 
around supermassive black holes in Active Galactic Nuclei 
(AGN). However, a detailed comparison with observations is 
limited by the lack of detailed models of self-gravitating discs 
and by an incomplete understanding of the basic physical 
processes involved. 

The importance of the disc self-gr avity in observation s 
of protostellar discs was pointed out bv lAdams et alj ([1988 ), 
who realized that a massive disc, with a flat rotation curve, 
could reproduce the observed flat long-wavelength Spectral 
Energy Distribution (SED) of many protostellar sources. 
However, in these early studies no detailed model of self- 
gravitating accretion discs was available, so that the required 
disc mass turned out to be unreasonably large. The se ideas 
were recently revisited by lLodato fc Bertin f <|200 lT) in th e 
light of a self-consistent disc model feertin fc Lodatoll9 99 ) , 
and they were able to model the observed SED of FU On- 
onis objects (a class of young stellar objects undergoing a 
phase of enhanced accretion), with a disc mass comparable 
to that of the central object. This model assumes that the 
outer disc (beyond a few au from the central star) is self- 
regulated. However, non-self-gravitating a-models of accre- 
tion discs would predict a rapidly declining radial profile 
of Q, which would eventually become unphysically small in 
the outer regions of the disc. In the case of FU Orionis discs 
the disc is predicted to be marginally s table already at a 
distance of ~ 1 au from the central star jBell fc Linlll994) . 
The arguments at the basis of the self-regulation mecha- 
nisms would suggest that in the outer regions of the disc 
an additional heating sour ce is required in order to keep 
the disc marginally stable. lLodato fc B ertin (2001) argued 
that the difficulty with a-models arises because these vis- 
cous models assume that energy dissipation is determined 
locally, whereas gravitational instabilities would naturally 
act in a global way, leading to a modification of the stan- 
dard estimates of the viscous dissipation power. 

Analogous cons iderations also ho l d in the context of 
AGN discs. Here, lLodato fc Bertin] (I2003T) have shown 
that the non-Keplerian rotation curve t raced by water 
masers in the Seyfert galaxy NGC 1068 iGreenhill et alJ 
19961: iGreenhill fc Gwinnlll997l) can be reproduc ed by self- 
following 
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Lodato fc BertirJ l 200ll) , have mod- 



regulated disc models 
the same arguments of 
elled the long-wavelength SED of AGN discs, based on the 
requirement that Q as 1. The latter authors also recog- 
ni se the need for some a dditional heating, but, contrary 
to lLodato fc BertirJ (j200lT l . attribute it to some "external" 
source, namely nuclear fusion in stars embedded in the disc 
(see also lCollin fc Zahnlll999]) . 

Actually, the issue of locality of transport in self- 
gravitating ac c retion discs still remains an open question. 
iLin fc Pringld (Il987f) have suggested that the transport in- 
duced by self-gravity could be described within a viscous 
framework, and introduced a modified a-prescription, where 
a ~ Q~ 2 . In this way, a self-regulated disc would be char- 
acterised by a rat her large effective viscosity, with a « 1. 
On the other hand. lBalbus fc Papaloizoul l)l999}l have shown 
that the energy equation for self-gravitating discs cannot be 



put in the form of a diffusion equation, as required in a 
viscous scenario, and that the energy flux contains some ex- 
tra terms, associated with global wave transport, that are 
"anomalous" from the point of view of v i scosity . Similar 
ideas had also been suggested bv lShu et all Jl99ot) . 

In t his context, numerical simulati ons play a cen- 
tral role. ILauehlin fc Bodenheimerl l)l994T> have performed 
global, smoothed particle hydrodynamics (SPH) simulations 
of massive isothermal discs, concluding that the density evo- 
lution o f the disc could be reproduce d in terms of simple a- 
models. ILauehlin fc Rozvczkal (Il996l) . using 2D grid based 
simulations, have shown that in order to reproduce the den- 
sity evolution induced by gravitational instabilities an a co- 
efficient dependent on radius was needed. However, these 
simulations did not include a detailed treatment of the heat- 
ing and cooling processes in the disc, which have been shown 
to play a fund amental role in determining the outco me of 
the instability dPickett et al.ll200d: iNelson et afll200Cl) . 

iGammiel l)200l|) has performed local, shearing-sheet, 
zero-thickness simulations of self-gravitating discs, includ- 
ing a simple cooling term, and has concluded that a local 
description is adequate in such "razor-thin" discs. However, 
Gammie's simulations are not appropriate to test global ef- 
fects, since locality is set up from the beginning, and they 
are only valid for infinitesimally thin discs, while the typical 
distance over which gravitational instabiliti es are expect ed 
to travel scales with the disc thickness. iRice et alJ <l2003al lbl) 
have already shown using global, three-dimensional simula- 
tions, how global effects can be important in the dynamics 
of self-gravitating discs, in relation to the related issue of 
the fragmentation of a massive disc. 

In this paper we present the results of global, three- 
dimensional, SPH simulations of massive, cooling, non- 
magnetised discs. Our main purpose is to quantitatively 
determine the dissipation power provided by gravitational 
instabilities and to compare the results with the expecta- 
tions based on a viscous theory of discs, in order to assess 
the extent to which the transport induced by gravitational 
instabilities can be regarded as a local process. 

The paper is organised as follows. In Section|5|we sum- 
marise the basic transport properties of viscous and of self- 
gravitating discs, introducing the basic physical quantities 
involved. In Section[3]we describe the numerical setup of our 
simulations. In Section^Jwe show the results of our compu- 
tations. In Section |S] we discuss our results in comparison 
with previous investigations and in relation to observed sys- 
tems. In SectionHJwe draw our conclusions. In the Appendix 
we discuss the more technical issue of transport induced by 
artificial viscosity in the numerical simulations presented. 



2 TRANSPORT IN ACCRETION DISCS 

In this Section we will summarise a few well-known results 
about the dynamics of viscous and self-gravitating accretion 
discs, that are going to be essential in the description of our 
results. We will not go into the full details of the derivation of 
the main result s, for wh i ch one can refer to standard reviews, 
such as that of lPringld Jl98ll) . 
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2.1 Non-self-gravitating discs 



2.2 Self-gravitating discs 



In the analysis performed in this paper we will adopt the 
thin disc approximation, and will therefore deal with verti- 
cally integrated quantities. The equations of motion for an 
axisymmetric disc, in cylindrical coordinates, are the equa- 
tion of continuity: 



9E 1 d , 



0, 



(2) 



and the azimuthal component of Euler's equation (expressed 
in terms of angular momentum conservation): 

where u is the radial velocity and Trj, is the relevant com- 
ponent of the viscous stress tensor, integrated in the verti- 
cal direction. This last term is the basic ingredient in the 
theory of accretion discs. As we have already anticipated, 
standard hydrodynamical viscosity is not sufficient to pro- 
vide accretion at the required rates, and therefore Tr^ is 
generally descr ibed by means of ad hoc prescriptions. The 
a-prescription JShakura fc Sunvaevlll97fl . based on simple 
arguments on the relevant physical scales of turbulent cells 
in the discs, assumes that is proportional to the disc 
pressure: 



dlnQ, 



din J? 



aEc» 



(4) 



where the proportionality constant a is an unknown pa- 
rameter, usually considered to be smaller than unity. The 
a-prescription can be also put in the following equivalent 
form, which involves the kinematical viscosity coefficient v: 



v — ac s _ff, 



(5) 



where H is the thickness of the disc. 

The effect of viscosity on the energy balance is twofold: 
viscous torques convect energy between neighbouring annuli 
of the disc and they dissipate energy. The energy which is 
convected across a ring at radius R per unit time is given 
by: 

dE 



dt 



while the dissipation rate per unit surface is given by: 
D(R) =T R4 ,\RU'\. 



(6) 



(7) 



If the disc is in thermal equilibrium, we can derive a useful 
relation between the viscosity coefficient a and the cooling 
timescale. If we assume that cooling can be simply parame- 
terised in the following way: 



£ C 2 



V 



icool 7(7 — l)*cool ' 



(8) 



where U is the internal energy per unit surface and 7 is 
the ratio of the specific heats, then, equating the viscous 
dissipation term, expressed in Eq. (J7J to the cooling term, 
expressed in Eq. JS), leads to: 



dlnfi 



dlni? 



7(7 - l)t cool Q' 



(9) 



where we have used also Eq. P). 



We will now turn to the transport properties of self- 
gravitating discs. The gravitational potential due to the disc 
will be denoted by & s , and g = — V$ s is the gravitational 
field. 

It can be shown llLvnden- Bell fc Kamaill972D that the 
equation of angular momentum conservation can be put in 
a form analogous to Eq. (JSJ, as required in a viscous frame- 
work, where the gravitational stress tensor is: 



1 R<t> — 



dz 



9R94> 

4ttG' 



(10) 



Eq. 11LH only accounts for the transport induced by the grav- 
itational field itself. Gravitational instabilities will also pro- 
duce density and velocity perturbations that contribute to 
the transport and should be included in the calculations. 
This contribution (the "Reynolds" stress) can be expressed 



(11) 



where 5v = v — u is the velocity fluctuation, with v the fluid 
velocity and u the mean fluid velocity. 

Given an expression for the gravitational stress tensor, 
one could therefore be tempted to use the a-prescription, 
and to assume that Tr^ is simply proportional to the local 
pressure. If large-scale structure and global processes play a 
role in self-regulating the disc, it could even be possible to 
give a "generalised" a-prescription, where a is to be de- 
termined by some global requirement (on this point, see 
lBertinlFl997l : ICoppil ll98C0. In any case, even if such global 
parameterisation is possible, it should be emphasised that 
the previous comments only apply to angular momentum 
transport. We still have to face the issue of energy trans- 
port in self-gravitating discs. In particular, we should check 
whether gravitational instabilities transport energy between 
neighbouring annuli according to Eq. © and whether they 
diss ipate energy according to Eq . £J. 

iBalbus fc Papaloizoul lll999T) have shown that in gen- 
eral energy transport cannot be described viscously. The 
energy balance equation for self-gravitatin g discs contains 
some "extra" terms (see Eq. (59) in IBalbus fc PapaloizovJ 
1999), that Balbus & Papaloizou ascribe to global wave 
transport. The important issue at this stage is to determine 
how important these extra terms are, and under which con- 
ditions they are able to influence the energy balance in the 
disc. 

In this work we will use global numerical simulations in 
order to compute explicitly the different physical quantities 
described in this Section. In particular we will evaluate the 
gravitational stress tensor, according to Eqs. 11UI and 1111 . 
and the corresponding "viscous" dissipation term, which can 
be then directly compared to the power actually dissipated 
in our simulated discs. 



3 NUMERICAL SETUP 
3.1 The code 

The three-dimensional simulations presented here were per- 
formed using smoothed particle hy drodynamics (SPH), a 
Lagrangian hydrodynamic code (see lBenj|l99ol : iMonaehanl 
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1992). In these simulations the central object is modelled 
as a point mass onto which gas particles can accrete if they 
approach to within the sink radius, while the gaseous disc is 
simulated using 250,000 SPH particles. Both the point mass 
and the gas particles use a tree to determine n eighbours and 
to calculate gravitational forces feend 1990'). and the cen- 
tral object is free to move under the influence of the disc 
gas. A significant saving in computa tional time is ma de by 
using individual, particle time-steps feate et al.l ll995') with 
the time-steps for each particle l imited by the Co urant con- 
dition and by a force condition (lMonaghanlll992l) . 

Since the main aim of this work is to check the en- 
ergy processes associated with gravitational instabilities, we 
explicitly solve the energy balance equation for the gas. We 
allow the disc to heat up due to both PdV work and viscous 
dissipation. The ratio of the specif ic heats is assum ed to be 
7 = 5/3. For the cooling, we follow Gammic (2001) and add 
a simple cooling term to the energy equation. Specifically, 
we use the same kind of parametrisation as in Eq. JHJ. For 
a particle with internal energy per unit mass u% , the cooling 
is implemented using 



with t cool = /3n~ MGammiel feOOlb and lRice etafl J2003bh 
have shown that for small cooling times the disc may frag- 
ment into gravitationally bound objects while for longer 
cooling times, the disc settles into a quasi-steady state 
with the imposed cooling balanced by dissipation through 
the growth of the gravitational instability. In particular, 
iGammiel lUpOl.) has shown th at the fragmentation boundary 
occurs for t coo i 3f2 _1 , while iRice et alj (|2003bj) show that 
global effects may lead to an enhanced tendency for frag- 
mentation, so that the critical value for t coo \ is somewhat 
increased. The reason why a disc that cools too rapidly is 
more prone to fragmentation is because gravitational heat- 
ing of the disc occurs on the dynamical timescale, so if the 
disc cools very rapidly, self-gravity has not enough time to 
prevent the formation of bound objects in the disc. In our 
simulations we have adopted (3 = 7.5, in which case none of 
our discs should or was found to fragment. 

3.2 Initial conditions 

We consider a system comprising a central star, modelled as 
a point mass with mass M 4 , surrounded by a disc with mass 
Miisc. We have performed three simulations with mass ratios 
q = M disc /M*, of 0.05, 0.1, and 0.25, respectively. The initial 
surface density profile is taken to be a power-law E oc R" 1 , 
while the initial temperature profile was T oc R~ 1/2 . The ini- 
tial velocity profile was calculated by including the enclosed 
cylindrical mass when determining the angular frequency Q. 
With these initial conditions, the minimum value of Q is at- 
tained at the outer edge of the disc. For each simulation the 
temperature normalisation is chosen such that the minimum 
value of Q is Q m i n = 2, and the whole disc is initially gravi- 
tationally stable. The disc is assumed to be in vertic al hydro- 
static equilibrium (see, for example. iPringle 1981), and we 
compute the scale height, H , using H — c s /Q and distribute 
the particles such that the vertical density profile is a Gaus- 
sian. Actually, in a self-gravitatin g disc, the vertical den sity 
profile is not rigorously Gaussian tBertin fc Lodatoll999h , so 



our initial setup, strictly speaking, is not in dynamical equi- 
librium. However, dynamical equilibrium is achieved rapidly 
(i.e., in a dynamical timescale) during the simulation. 

Our calculations are essentially scale-free. In dimension- 
less units, the disc extends from R[ n = 0.25 to -Rout = 25, 
and we have taken M* = 1. In these units, one dynamical 
timescale (orbital period) at R = 1 is equal to 2ir code units. 
Therefore, one orbital period at the outer edge of the disc is 
roughly equal to 800 time units. 

The disc is not initially in thermal equilibrium since 
the main source of heating, i.e. gravitational instabilities, is 
turned off because the disc is initially stable, while cooling 
is already effective. We follow our simulations for ~ 5000 
time units (i.e. approximately 6 orbital periods at the outer 
edge of the disc) , by which stage the whole disc has reached 
thermal equilibrium. 

3.3 Artificial viscosity: the Balsara switch 

The standard SPH viscosity (eg.. iMona ghan 1992) consists 
of a quadratic term similar to a Von Neumann-Richtmeyer 
viscosity (characterised by a dimensionless coefficient called 
/3sph), and a linear term that introduces a viscosity in shear- 
ing flows (characterised by a dimensionless coefficient called 
Qsph). Since the goal of this work is to study transport 
associated with gravitational instabilities, we would like to 
reduce any angular momentum transport associated purely 
with the artificial viscosity (i.e., we wish to reduce the arti- 
ficia l shear v i scosity ) . 

iBalsaral lll995l) suggested adding a shear correction 
term, known as the Balsara switch, to the standard SPH 
artificial viscosity which maintains the viscosity in compres- 
sional flows but reduces it to zero in pure shear flows. For 
all the simulations presented here we have used the Bal- 
sara form of the artificial viscosity and have used a value of 
0.1 for the coefficient of the linear artificial viscosity term 
(qsph) and a value of 0.2 for the coefficient of the quadratic 
viscosity term. A detailed discussion of the transport of an- 
gular momentum induced by the chosen artificial viscosity 
is presented in the Appendix. Here we anticipate that the 
angular momentum transport induced by artificial viscosity 
is at least a factor 10 smaller than that induced by grav- 
itational instabilities in our simulated discs, and therefore 
plays a minor role in the dynamics of the disc. 



4 RESULTS 

4.1 General features 

In all our simulations the disc initially cools down until the 
value of Q becomes small enough for gravitational instabil- 
ities to become effective and to provide a source of effective 
heating. At later stages the disc settles into a quasi-steady 
state, where the heating provided by the instabilities bal- 
ances the imposed cooling term. As predicted by the ar- 
gument of self-regulation, this quasi-steady state is charac- 
terised by an almost constant value of Q ~ 1, over a wide 
radial range. The main features of the simulations are sum- 
marised in Figs. El 013 andgl 

Fig. shows the surface density structure for the three 
simulations, with different total disc mass (q — 0.05, q = 0.1 
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and q — 0.25). The global structure induced by self-gravity 
is clearly seen in all three cases. The density enhancement 
in the spiral arms Ap/p typically ranges between 2 and 4 
in all cases. Already from this figure it can be noticed that, 
as the total disc mass increases, the pattern of the instabil- 
ity becomes progressively more dominated by low-m modes. 
This is confirmed by the Fourier analysis shown in Fig. [5] 
that shows how in the higher total mass cases the modes 
with m < 5 dominate the structure. The Fourier amplitudes 
in Fig. [5] are computed as follows: we divided the disc in 
concentric rings of width A_R = 2.5 in dimensionless units. 
For each ring we then compute the mode amplitude A m as: 



E 



N ring 



(13) 



where N r i ng is the number of particles in the ring, and </>; is 
the azimuthal angle of the z-th particle. 

Fig. shows the radial profile of Q at three different 
times, towards the end of the simulation. The profile is not 
significantly altered with time, indicating that the simula- 
tions have reached thermal equilibrium, and are in a quasi- 
steady state. Actually, the fact that the self-regulated value 
of Q turns out to be so close to unity in our simulations 
is quite remarkable. In fact, the marginal stability value Q 
is equal to unity only when the perturbation analysis is re- 
stricted to the case of the axisymmetric stability of an in- 
finitesimally thin disc. The disc thickness has an important 
stabilising effect, which leads to a lower temperature (and 
thus a lower value of Q) of the marginally stable state; on 
the other hand, even a light disc can be destabilised by the 
presence of swing amplified non-axisymmetric disturbances 
with high m (for a discussion of these mechanisms see, e.g., 
lBertinll2000h . It then appears that these two competing ef- 
fects basically counterbalance each other. A detailed inves- 
tigation of the combined effect of these physical mechanisms 
in determining the precise value of Q at marginal stability 
is, however, beyond the scope of the present paper. 

Fig. [I] displays the initial and final radial profiles of the 
surface density (averaged in the azimuthal direction). It can 
be noticed that there is no significant evolution of the sur- 
face density at large radii. The major changes in the surface 
density occur close to the boundaries. In particular, there is 
a rapid drop in surface density close to the inner boundary. 
This is due to the fact that we have not attempted to give a 
detailed description of the boundary layer. The sudden lack 
of SPH particles at R < 0.25 causes an increased artificial 
pressure which pushes the inner particles into the sink radius 
of the star. This effect should be reduced with a more accu- 
rate description of the boundary layer. We have also checked 
that the total angular momentum of the disc is conserved to 
within a few percent throughout all the simulations. 

When low values of the artificial viscosity are used, par- 
ticle interpenetration might lead to a poor representation of 
strong shocks in SPH. This is not a serious issue in our case, 
since in our simulations only mildly supersonic shocks are 
involved. Based on the density contrast in the spiral arms, 
we estimate the Mach number of the shocks to be M < 1.5. 
At these low values of M, a value of Psp h ~ 0-2 is a lready 
sufficient to stop particle interpenetration jBatell995D . This 
is confirmed by the well defined spiral structure that we 
obtain (that would have been smeared out if significant par- 



ticle interpenetration was indeed present), consistent with 
the results of previous simulations that used the standard 
SPH viscosity and higher values for the viscosity coefficients 
jRice et al]l2003allrl . 



4.2 Angular momentum transport and energy 
dissipation 

The torque produced by gravitational instabilities in the disc 
is given by the sum of the two terms described in Eqs. 1101 
and 111 It . After averaging the stress tensor azimuthally and 
radially, over a small region AR = 0.1, we compute the 
corresponding value of a (see Eq. J1J): 



ct(R) 



dlnfi 



dlnR 



1 < Tgr > + < j£; yn > 

Sri 



(14) 



The resulting radial profiles of a are shown in Fig.|5|for 
the three cases q = 0.05, q = 0.1, and q = 0.25. The upper 
panels show separately the hydrodynamic and gravitational 
contributions to a, while the bottom panels show the sum 
of the two. The plots show the time average of a at the end 
of the simulation, once the disc has reached a quasi-steady 
state. The time-averaging interval is 500 time units, i.e. 0.6 
orbital times at the outer disc edge. 

We can now use the general results of viscous disc the- 
ory outlined in Sec. [5] to test the locality of transport in 
our simulations. In fact, Eq. © gives us firm expectations 
for the value of a needed to balance the imposed cool- 
ing, if energy dissipation can be treated in a viscous frame- 
work, i.e. by using Eq. 0. In particular, in our simulations 
i coo ifi = /3 = 7.5, 7 = 5/3 and, since our discs are nearly 
Keplerian, dlnfi/dlni? ~ —3/2. Inserting these numbers in 
Eq. JUJ would give us an expected value of a ~ 0.05. It is 
important to note that the fact that the expected a turns 
out to be nearly independent on radius is a result of choosing 
the cooling time to be simply proportional to the dynami- 
cal time-scale. In the general case, of course, the resulting a 
needs not be constant. The dotted line in Fig. shows the 
expected value of a. Our results are in fairly good agreement 
with the expectations of viscous transport theory. 

It is also interesting to compare the dissipation rate 
D(R) that would result if the transport process were vis- 
cous, i.e. computing D(R) based on Eq. 0, with the actual 
dissipation rate. Since our discs are in thermal equilibrium, 
the dissipation rate can be obtained from the cooling rate 
using Eq. @. This comparison is shown in Fig. HJ where 
the solid line shows D(R) computed from Eq. 0, while the 
dotted line shows the cooling rate. Again, the plots show a 
good agreement with the expectations, the only exception 
being the outer parts of the most massive disc, where the 
actual dissipation rate seems to be slightly larger than that 
expected from a viscous process. 

In general, our results indicate that, up to disc masses 
of Afdisc = 0.25M*, gravitationally induced transport is rea- 
sonably well described within a local framework. This result 
could also be anticipated since the disc dynamics in all the 
three simulations are dominated by rather high-m modes, 
that dissipate on a short length-scale. 

As a separate test for the locality of the transport, we 
have also computed a part (-R, d), which we define as the grav- 
itational part of a(R), taking into account only those parti- 
cles inside a spherical radius d from the radial point R where 
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we compute the stress. This quantity gives us a measure of 
the size of the region that has the most significant contri- 
bution to the gravitational stress at a given point. However, 
it should be kept in mind that this quantity only accounts 
for the stress produced by the gravitational field, without 
including the hydrodynamical component. 

Fig. Qshows the results for R — 15 for the three simula- 
tions. The upper panel shows clearly that the more massive 
the disc, the more distant regions from the point contribute 
to the stress. However, it should be noted that, since in the 
thermal equilibrium state Q » 1 in all the three 
more massive disc is also characterised by a larger value of 
the effective thermal speed c 3 , and is therefore thicker. Since 
we expect that gravitational disturbances propagate over a 
length-scale of the order of the thickness of the disc, to com- 
pare appropriately the results of the three simulations, we 
should check the dependence of a par t on d/H. This is shown 
in the bottom panel of Fig. [7| In all cases, more than 80% 
of the contribution to the gravitational torque comes from 
a region with size A « 10H. Transport is local if A <C R. 
We can therefore conclude that angular momentum trans- 



port is determined by local conditions only for discs with 
H/R <C 0.1. When Q ~ 1, the ratio H/R is proportional to 
Mdisc/M*, which means that sufficiently massive discs may 
violate the previous condition. Actually, our most massive 
disc (for which M,n BC /M+ = 0.25) only marginally satisfies 
H/R < 0.1. Indeed, as can be seen from Figs.|^|and|S| global 
wave transport might play some role in the outer disc, in this 
case. 



5 DISCUSSION 

5.1 Comparison with previous work 

The distinctive feature of the present work is that we have 
performed global, three-dimensional simulations of massive 
discs, including the detailed effects of heating from gravita- 
tional instabilities and cooling. Ou r simula tions are similar 
to thos e performed by iRice et all i2003al) and iRice et all 
(2003b), with the difference that these previous investiga- 
tions were concerned about either the motion of the central 
object caused by the massive disc, or the issue of fragmen- 
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Figure 2. Amplitude of the first Fourier components of the density structure for different radial ranges in the three cases: (upper left) 
q = 0.05, (upper right) q = 0.1, (bottom) q = 0.25. 



tation of the disc, while here the main goal is to characterise 
the transport properties induced by gravitational instabili- 
ties. 

We believe that both performing 3D simulations and 
solving explicitly the energy equation are essential to deter- 
mine the final outcome of the instabilities. In fact, the dy- 
namical properties of self-gravitating discs are determined 
to a large extent by the process of self-regulation, which 
is strongly dependent on the detailed heating and cooling 
mechanisms, as outlined in the Introduction. Furthermore, 
since one of the main tests we want to perform is to check 
the kind of dissipative process associated with gravitational 
instabilities, solving the energy equation is essential. The 
requirement of 3D simulation is also fundamental, since the 
typical size of gravitational disturbances is related to the 
disc thickness, so that any zero-thickness simulation will au- 
tomatically lead to an underestimate of global effects. 

Previous numerical work done on the subject includes 
global, 3D SPH simulations o f massive isothermal discs 
jLauehlin fc Bodenheimerlll994l) ; global, 2D, SPH simula- 



tions with detailed heating and cooling (Ne lson et aljfe oOO): 
global, 3D, grid-bas ed simulations with heating and cooling 
JPickett et all [2000 ) ; and l ocal, 2D grid-b ased simulations 
with heating and cooling jGammid 1200 lTl . In this Section 
we describe the similarities and the differences between our 
study and these previous ones. 

One of the first stud ies of gravitational instab il ities i n 
discs was performed by lLaughlin fc BodenheimeU lll994) . 
They modelled a very massive disc (with Mdi sc ~ M*) and 
followed its evolution with a 3D SPH code, without includ- 
ing any heating or cooling term, but simply assuming that 
the disc was locally isothermal. In this study, the authors 
also tried to give a detailed characterisation of the trans- 
port. Their approach was however slightly different to ours, 
in that they evolved their simulation long enough to capture 
the viscous evolution of the disc, and then compared the 
evolution of the azimuthally averaged surface density with 
simple one-dimensional viscous models, concluding that the 
disc evolution could be well reproduced by a viscous model 
with a m 0.03. This work is important because it clarifies 
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Figure 3. Profiles of the Q parameter for the three simulations: (upper left) q = 0.05, (upper right) q = 0.1, (bottom) q = 0.25. 



that gravitational instabilities are actually able to transport 
angular momentum efficiently and that the surface density 
evolution is indeed of a diffusive nature, as expected (see 
Section but does not answer the important question of 
whether energy dissipation is local or global. 



iNelson et alJ i2000l) performed 2D simulations with par- 
ticular emphasis on the cooling processes in the disc and 
included a more realistic cooling function than the sim- 
ple parameterisation adopted here. Their disc mass was 
Afdisc = 0.2M*, very similar to our most massive case. They 
estimated the effective a associated with artificial viscosity 
(see Section I3.3H to be of the order of 5 x 10~ 3 , compara- 
ble to the expectations based on Eqs. IA1I and 1A2L and 
an order of magnitude smaller than what we get here from 
transport induced by gravitational torques. Indeed, in our 
simulations most of the angular momentum is carried by col- 
lective instabilities and by gravitational torques, rather than 
by the artificial shear. Comparing the dissipation rates from 
shock he ating with that from turbulent heating in their sim- 
ulations, Nelso n et al.l i2000f) conclude that, at least in the 
outer disc, gravitational torques should not contribute more 



than the torques associated with artificial viscosity, while in 
our simulation gravitational torques are everywhere domi- 
nant. A possible reason for this discrepancy lies in the 3D 
nature of our simulations, which may allow gravitational 
disturbances to travel further and allows low-m modes to 
be more prominent. Actually, even if the mode amplitudes 
that we obtain her e (see F ig. El arc in rough agreement with 
what obtained bv lNelson et alJ i200(J) . they find very sim- 
ilar mode amplitudes for all modes with m < 8, while in 
our case there is a marked increase of mode amplitude for 
modes with m < 5 (see Fig. 

iGammiel |2001) performed an analysis very similar to 
our own (actually, we adopt the same prescription for the 
cooling term, see Eq. 1121 1. The main difference is that Gam- 
mie's simulations are local and two-dimensional while ours 
are global and three-dimensional. He computes the effective 
a based on Eq. 1141 . as we do. The main results that we 
obtain are basically in agreement with Gammie's results, in 
that the expected value of a from a viscous theory of discs 
is very close to the computed value in the simulations. How- 
ever, by using a global approach, we are now able to check 
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how these results depend on the total disc mass and on the 
disc thickness, whereas Gammie could only extrapolate his 
results to thi cker configurations. 

Re cently. I Johnson fc Gammid i2003tl have extended the 
work of iGammiel ^OOlF to include a more realistic cooling 
function in 2D simulations. Actually, the use of detailed cool- 
ing functions makes it even more important to perform a 
full 3D simulation, since most of the radiative transport will 
occur in the vertical direction, so that using vertically aver- 
aged values for the relevan t physical quantities, as done by 
Ijohnson fc Gammid <l2003f) . might lead to incorrect results. 

5.2 Impact on observed systems 

It is now commonly accepted that most low-mass young 
stellar objects possess a circumstellar disc, with lifetimes 
of at least lMyr. Circumstellar discs play an important 
role in the process of star formation and are believed 
to be the site where planet formation occurs. Observa- 



tional evidence for the presence of such discs is either in- 
direct, i.e. based on the disc emission at long- wavele ngths, 
such as in the infrared (star ting from [ A dams et alj Il988l) 
and at sub-mm wavelengths feeckwitlr^tHuT^jQOIK or di- 
rect, Le ; _l2v_jmaghig_of_siWwMe^ discs in the Orion neb- 
ula llMcCaughrean fc O'Delil I1996T) . Especially in the ear- 
liest phases of star fo rmation these discs m i ght b e fairly 
massive: for example, lLaunhardt fc Sargend (1200 ll) . using 
sub-mm observations, report the discovery of a massive 
disc (with Mdi sc /M* > 0.3) in a very young (Class I) 
protostellar object. There are also some indications that 
the discs in FU Orionis ob jects might be fairly massive: 
ISandell fc Weintraubl (|200^ report M d i ac /M* > 0.1 in most 
of the FU Orionis discs they have observed. As already men- 
tioned, detailed modelling of FU Orionis outbursts produce 
radial profiles of Q that fall below unity already at a dis - 
tance of ~ 1 au from the central object ijBell fc Linlll994F ). 
All these systems are likely to be affected by self-gravity 
that, if indeed energy dissipation is non-local, may produce 
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Figure 5. Effective a produced by gravitational instabilities for (upper left) q = 0.05, (upper right) q = 0.1, and (bottom) q = 0.25. 
The top panel shows the separate contribution of a gra v and OR eyn , the lower panel shows the sum of the two contributions compared 
with the expected value from a local viscous model (dotted line). 



some observable modification in the SEP. It is th en inter- 
esting to see that the models bv lBell fc Linl dl994l) predict 
H/R > 0.1, which, according to the results of this work, are 
large enough for non-local effects to become important, as 
suggested bv iLodato fc Bertinl J200ll) . 

In the context of AGN discs, the distance at which the 
disc becomes marginally stable to gravitational instabilitie s 
is typically of the order of 10 3 7? g iLodato fc Bertin| [2003). 
where R s is the g ravitational radius of the b lack-hole. Water 
maser emission iGreenhill fc Gwin nlll997l) and radio con- 
tinuum observations llGalllmoTe^^ilT^99^ show that in 
many cases the disc can extend to radii much larger than 
that, thus allowing self-gravity to influence the disc struc- 
ture. In this context, self-regulated models have been ap plied 
both to the modelling of the SED dSirko fc Goodmanll2003l) 
and to the modelling of the rotatio n curve in the outer dis c 
of the Seyfert galaxy NG C 1068 (ILodato fc Bertinl l2003t) . 
I.Iohnson fc Gammiel l|200.?l) use their 2D results with "real- 
istic" cooling (see previ ous Subsection) t o argu e that the 
disc model proposed by ILodato fc Bertinl teOO&t) for NGC 



1068 would be subject to f ragmentation, but unfortunately 
Ijohnson fc Gammiel i2005t) do not explore the region of the 
parameter space relevant to the Lodato & Bertin model. 

As a final comment, we note that the present work refers 
to the situation where the dominant source of transport and 
dissipation is provided by gravitational instabilities. In ob- 
served systems, other sources of transport could be present, 
which might reduce the effect of gravitational instabilities. 



6 CONCLUSIONS 

In this paper we have investigated the transport properties 
induced by disc self-gravity in relatively massive accretion 
discs. In particular, we have discussed the extent to which 
angular momentum transport and energy dissipation can be 
described within a viscous framework. To this goal, we have 
performed global, three-dimensional simulations using SPH. 

Most of the recent work on the subject has focused on 
the issue of the conditions for disc fragmentation, especially 
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Figure 6. Viscous dissipation rate (solid line, from Eq. J7J and actual cooling rate (dotted line) for (upper left) q = 0.05, (upper right) 
q = 0.1, and (bottom) q = 0.25. 



in relatio n to the process of planet formation in protoste l- 
lar discs <Bosdl2002l: iRice et alJl2003bt iMaver et al.ll2002t) 

or massive star fo rmation in AGN discs jGoodman fc Tar] 
2003i; iLevinl 120031) . Here, we would like to check the be- 
haviour of self-gravitating discs in the case where fragmen- 
tation does not occur and the disc evolves towards a quasi- 
steady state where cooling is balanced by heating from grav- 
itational instabilities. 

The issue of the transport properties induced by grav- 
itational _instabifitiej3_Ju^^ analytically in the 
pa st by ll^yrid^n^BelLfc^ajiiaisj lll972l) and, more recently, 
bv lBalbu^^Papaforzouri^99^ . There are two major as- 
pects of the problem: (i) to what extent is the shear stress 
Tb.4, at a given radius R determined by local conditions? This 
point directly calls into question the use of the a formalism, 
which explicitly requires that Tr,^ is only dependent on the 
local values of density and thermal velocity (see also the 
discussion in Section r2.21 . The second aspect is (ii) whether 
energy transport is only determined by the gravitationally 
induced shear stress or whether instead global wave trans- 



port occurs, hence influencing energy diss ipation in the disc, 
as argued bv lBalbus fc Papaloizoul (Il999l) . 

We have performed several simulations, with different 
values of the ratio M^isc/M*, to test how the global trans- 
port properties of the disc change with increasing total disc 
mass. Our simulations appear to evolve to a quasi-steady 
state, characterised by an almost flat profile of the axisym- 
metric stability parameter Q « 1, in which the heating 
provided by gravitational instabilities balances the imposed 
cooling rate. We have computed the torque induced by grav- 
itational instabilities and the corresponding viscous dissi- 
pation rate, which we then compare to the actual dissipa- 
tion rate in our simulated discs. We have found values of 
a ~ 0.05, in reasonable agreement with the expectation from 
a viscous theory. Energy dissipation in our simulations is also 
fairly well described using a viscous approach. These results 
directly address aspect (ii), described above , and confirm 
the argument of lBalbus fc Papaloizoul il999l) . who had in- 
deed argued that for "self-regulated" discs, in which Q ~ 1, 
global wave transport of energy would not play a major role. 
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On the other hand, we have indeed noticed that more 
massive discs tend to be dominated by lower m modes, 
leading to a more global pattern of the gravitational dis- 
turbances. We have also been able to determine the size of 
the region that mostly contributes to the torque at a given 
point in the disc. More than 80% of the torque is produced 
within a region of size A ~ 10H, where H is the thickness of 
the disc. These results directly address issue (i), described 
above. We can therefore conclude that a viscous description 
of the transport in self-gravitating discs is only appropriate 
when H/R < 0.1. In system s like FU Orionis , where models 
predict disc thickness ~ 0.1 (Bell fe Lin 1994), global effects 
could play a role and modify the dissipation rates in the 
outer disc. 

The results of this work should be extended both to- 
ward a more thorough investigation of the parameter space 
(in particular, it is very important to test the transport prop- 
erties of discs more massive than 0.25M*, and to explore the 
results obtained with different initial conditions and cooling 
timescale, which in the present work has been taken to be 
icooi = 7.5r2 _1 ), and toward a more real istic description of 
the cooling function (see, for example, iNelson et alJ l200ol 
and Ijohnson fc Gammiell200of) . 
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Figure Al. Radial surface density plots at times of t = 0.0, t = 0.4, and t = 1.0 for an axisymmetric ring evolved using SPH with the 
standard artificial viscosity (dots) and evolved analytically (solid lines) with a viscosity determined using the parameters of the SPH 
calculation. 




Figure A2. Radial surface density plots at times of t = 0.0, t = 0.4, and t = 1.0 for an axisymmetric ring evolved using SPH (dots) 
and evolved analytically (solid lines). In the left hand figure the Balsara switch was used in the SPH artificial viscosity while in the right 
hand figure the simulation was evolved using the standard SPH viscosity. The Balsara switch clearly reduces the effective shear viscosity 
by a factor of between 5 and 10. 



APPENDIX A: ANGULAR MOMENTUM 
TRANSPORT INDUCED BY ARTIFICIAL 
VISCOSITY 

In this Appendix we discuss the magnitude of the angular 
momentum transport associated with the artificial SPH vis- 
cosity, in order to be sure that the main contribution to the 
shear stress described in Section 2] is actually due to the 
effect of gravitational disturbances. 

Al The Balsara switch 

As mentioned in the main text, we have adopted the Balsara 
form of viscosity in order to reduce the effect of artificial 
viscosity in transporting angular momentum. Although this 
modification to t he standard SPH viscosity has been studied 
in so me detail dNavarro fc Steinmet d IT997t iThacker et alJ 
2000), it would be useful to have some idea of the reduction 



in shear viscosity that occurs when using this form of the 
artificial viscosity. 

In the continuum limit, the linear term in the standard 
SPH artificial viscosity can be shown to have th e following 
form jArtvmowicz fc LubowllT99llMurravlll99ri) 

v = iaspHCs/i (Al) 

o 

where cksph is the linear viscosity coefficient, c s is the sound 
speed and h is the SPH smoothing length (which essentially 
determines the resolution of the simulation). As shown by 
iLvnden-Bell fc Pringld lll974l) , the time evolution of the sur- 
face density and radial velocity of an initial Gaussian ring 
can be d etermined anal ytically. To test the standard SPH 
viscosity, iMurravl dl996t) therefore considered the evolution 
of an initial Gaussian ring with the initial radial velocity 
determined analytically, using Eq. lAH to determine the 
viscosity. We have repeated this calculation and show the re- 
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suits in Figure lATl The SPH calculation had a linear viscos- 
ity coefficient of cisph = 10, a sound speed of c a = 0.02, and 
a constant smoothing length of h = 0.0075. The choice of 
the SPH param eter wa s made in order to be consistent with 
IMurravl jl99d) . As in iMurravl §996), the pressure forces 
were switched off so as to study the artificial viscosity in iso- 
lation. This, on the one hand, prevents the ring from spread- 
ing due to pressure forces, but on the other hand would 
cause the disc to collapse to the midplane. As a check that 
converging flows do not influence our results, we have also 
performed some tests with a fixed vertical coordinate of the 
SPH particles, and we found no significant differences. The 
dots in Figure lATl show the time evolution of the SPH sur- 
face density. For the chosen SPH parameters, the associated 
shear viscosity, according to Eq. 1A1I . is v = 1.9 xlO" 4 . The 
solid lines in Figure |A"T1 show the analytic viscous evolution 
of an initial Gaussian ri ng with an imp osed shear viscosity 
of v = 1.9 X 10~ 4 . As in IMurravl Jl996l) the SPH evolution 
of the initial Gaussian ring follows the analytic result very 
closely and Eq. HAH appears to be a good representation 
of the shear viscosity associated with the linear term in the 
standard SPH artificial viscosity. 

To perform a similar calculation to determine the vis- 
cous transport when using the Balsara switch is more subtle 
since the shear viscosity should, ideally, go to zero. This 
will, of course, not be exactly true in practice, but we can- 
not now determine what the initial radial velocity profile 
should be. To get some idea of the viscous transport when 
using the Balsara switch we have considered the time evo- 
lution of a Gaussian ring in which, since ideally we expect 
no shear viscosity, we set the initial radial velocities to zero. 
For linear viscosity coefficients of qsph = 0.1 and lvsph = 1 
there is no noticeable spreading between t = and t = 1. 
For asm = 10, however, the ring did spread slightly. The 
left hand side of Figure IA2I shows the SPH evolution of an 
initial Gaussian ring (dots) plotted at times of, as in Fig- 
ure E3 t = 0.0, t = 0.4, and t = 1 for a S p H = 10. A 
direct comparison with the analytical prediction for the sur- 
face density evolution and with the results shown in Fig. lAll 
would be misleading in this case, because initially the ring 
has to settle down before attaining the appropriate radial 
velocities resulting from the viscous evolution. As a result, 
the initial e volution of the ring is slower th an predicted ana- 
lytically bv lLvnden-Bell fc Pringlel dl974h . In the left panel 
of Figure IA2I we also plot the analytic evolution of an ini- 
tial Gaussian ring (solid line) with an imposed viscosity of 
v = 8.0 x 10~ 6 , which appears to best describe the average 
evolution of the ring. In order to compare these results with 
those obtained from the use of the standard SPH viscosity, 
we performed an identical calculation (with the initial radial 
velocities set to zero) except using the standard SPH viscos- 
ity without the Balsara switch. This is shown in the right 
hand side of Figure lA"2l Again the dots represent the evolu- 
tion of the SPH surface density at times of t — 0.0, t = 0.4, 
and t = 1 and the solid line shows the analytic evolution 
of the surface density. Also in this case, as discussed above, 
the initial evolution is somewhat slower than expected. The 
average evolution between t = and t = 1 can be fitted 
with a viscosity v = 4.5 x 10~ 5 . These results, though not 
conclusive, indicate that the Balsara switch reduces the ef- 
fective shear viscosity by a factor of between 5 and 10. By 
using the Balsara switch we should therefore be able to re- 
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Figure A3. Angular momentum transport induced by artificial 
viscosity in our SPH simulations, as parametrised by an effective 

«art- 

duce the angular momentum transpo rt due to the artificia l 
shear viscosity without, according to lThacker et all |2000), 
significantly influencing the handling of shocks. 

A comparison between Eq. IA1I and Eq. JSJ readily 
shows that the artificial shear viscosity in SPH (not using the 
Balsara swith) can be described in terms of the a parametri- 
sation in the following way: 

1 h , . . 

Oart = -Q!SPH-— ■ (A2) 
8 tl 

Eq. 1A21 therefore shows that the magnitude of a ar t depends 
on how well the vertical structure of the disc is resolved. In 
most of our simulations we typically had H « 5/i, so that, 
if we had used the normal SPH viscosity, we would expect 
Q art ~ 2.5 x 10 -3 , for asra = 0.1 (the value used in all our 
simulations). The use of the Balsara switch enables us to 
further reduce this contribution to a ar t ~ 5 x 10 -4 . 

A2 Artificial transport in disc simulations 

As discussed in the previous Section, thanks to the Balsara 
switch, we expect the linear term in the viscosity to produce 
only a minor contribution to the shear stress. However, the 
quadratic term of the artificial viscosity (parameterised by 
/3sph) might still give some contribution to the angular mo- 
mentum transport. In addition, in our self-gravitating disc 
simulations, we have also added pressure force, which might 
as we ll contribute to the stress, as found also by IMurravl 
(1996). Moreover, we would also like to test this issue in a 
less idealised case, with respect to the simple spreading ring 
calculations described in the previous Section. For this pur- 
pose, we have performed some simulations taking the same 
surface density profile as in Section 01 (i.e. E oc R' 1 ) and 
a constant Q profile, with Q = 1, but in which we have 
switched off the main source of transport and dissipation 
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(i.e. the disc self-gravity) and the cooling. The artificial vis- 
cosity coefficients were «sph = 0.1 and /3sph = 0.2, as in 
the simulations presented in Section 0] We have computed 
the Reynolds stress according to Eq. 1111 . and obtained an 
effective a art from Eq. The results are shown in Fig. 
IA3I The magnitude of this artificial transport is never larger 
than a art ~ 5 10 -3 . Since our results (see Sectional indi- 
cate that the transport induced by gravitational torques can 
be parametrised with an effective a~5x 10~ 2 , we can be 
confident that the major contribution to the stress tensor 
comes from gravitational torques rather than from artificial 
viscosity. 



